Monte Carlo study of linear chain submonolayer structures. 
Application to Li/W(112) 
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The lattice gas model for adsorption of alkaline elements on W(112) surface is studied by Monte 
Carlo simulations. The model includes dipole-dipole interaction as well as long-range indirect 
interaction. The numerical results show that truncation of the indirect interaction even at 200A may 
change a phase diagram, i.e., new phases containing domain walls might occur. It is demonstrated 
that a defected phase can exist at high temperatures even if it is not stable at T=0. The phase 
diagram for Li/W(112) is constructed and long periodic chain structures (9x1), (6x1), (4x1), (3x1), 
and (2x1) are found to be stable at low temperatures. Role of thermal fluctuation is discussed by 
comparison of Monte Carlo results with mean field approximation results. 
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PACS numbers: 68.35.Rh, 05.10.Ln 



I. INTRODUCTION 



Chemisorption on metal surfaces has been studied 
for over 30 yearsiiSiiiiSiSiLSiS. One of interesting prob- 
lems concerns with structures and phase transitions 
in chemisorbed submonolayers. In particular alkaline, 
alkaline-earth, and rare-earth elements adsorbed on 
W(II2) and Mo(112) form linear chain structuresii^i. In 
case of low coverages (O < 0.5), the commensurate 
phases with chains normal to the troughs of the substrate 
are observed. The period of these structures changes 
from 2 to 9i. At higher coverages (O > 0.5), incommen- 
surate surface phases occur^-^. It has been suggested^'- 
that formation of chain structures might be due to in- 
direct interaction between adatoms via conduction elec- 
trons. Recently we proposed^" the lattice gas model 
to account for linear chain structures on W(II2) and 
Mo(112). The model includes dipole-dipole interaction 
as well as indirect oscillatory interaction. As the inter- 
actions are of long range, the model is difficult to study 
theoretically. Preliminary studies performed in molecular 
field approximation confirmed that competition of dipole- 
dipole interaction and indirect interaction leads to forma- 
tion of linear chain structures. Further extensive analy- 
sis of ground states in related effective one-dimensional 
model with infinite range of interactiongii pointed out 
that competition between these two interactions is cru- 
cial in formation of linear chain structures. In very re- 
cent investigation of phase transitions in Li/Mo(II2) and 
Sr/Mo(II2), H. Pfniir et al^i^ have shown that the 
lattice gas model could be useful in description of long 
periodic phases observed at low coverages. Results of 
their investigation support an idea that surface electronic 
state are responsible for the oscillatory indirect interac- 
tion between adatoms. 

It is intention of this paper to study phase diagrams 
by Monte Carlo simulations. This method was suc- 
cessfully applied to lattice gas models with short range 
interactionsi^ii^iiSiii. On the other hand, one can expect 



more difhcult problems in case of long range interactions, 
e.g., occurrence of complicated structuresiSiiSiSSiSi. The 
paper is organized as follows. In section 2 the lattice gas 
model is described. Two difficulties connected with ap- 
plication of Monte Carlo simulations to this model are 
discussed in section 3. In section 4 phase diagrams are 
presented and section 5 contains discussion and conclu- 
sions. 



II. THE MODEL 

In the previous paper— we introduced the lattice gas 
model of linear chain submonolayer structures on the 
(112) surface of W and Mo. Adsorption sites form the 
rectangular lattice (see Fig. with = and 
Gy — a\/3/2, where a denotes the lattice constant of the 
bcc elementary cell. The model Hamiltonian within the 
grand canonical ensemble, is defined as 



n 



(1) 



where rii is the occupation variable at each adsorption 
site i with rii = 1 if the ith site is occupied by an adatom 
and rii — 0, if not. The chemical potential is denoted 



by fi and Vi 



V{r1 



J) is the pairwise interaction 



consisting of electrostatic and indirect interactions 
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for y ^ 



(2) 

The ffi'st part of Eq. Q, 2(f/\r\^ , describes a repul- 
sive dipole-dipole interaction between two adatoms with 
identical dipole moments d. The second term in Eq. Q 
represents an asymptotic part of the indirect interaction 
between adatoms via conduction electrons of the sub- 
strate. This particular interaction is highly anisotropic 



2 



• R • ^ 


' E • \ • 

b \ 






m 







x=[ 1 1 ] 



FIG. 1: A lattice of adsorption sites on (112) surface of bcc 
metal. The elementary cell is drawn in the left corner. Ranges 
of two-particles interactions used in simulations are depicted. 
For a notation see the text. 



and is closely related to the existence of nearly flattened 
segments of the Fermi surface being perpendicular to 
the [ill] axis directed along the atomic troughs of the 
substrat o^'^? . It is also possible that contribution to 
this interaction comes from quasi-one-dimensional sur- 
faces statesH. Very recently, G. Godzik et al.^"^ suggested 
that oscillatory indirect interaction of this form is, in case 
of Sr/Mo(112), induced by charge density waves involv- 
ing surfaces states. An amplitude A and a phase if can 
be treated as phenomenological parameters and kp de- 
notes a wavevector of electrons at the Fermi surface. Fi- 
nally, for y = there is an attractive indirect interaction 
Eh between the nearest neighbours along the x-direction 
which is responsible for the chain formation. 

In this work we choose the set of model parameters 
the same as proposed in the previous paperi? for adsorp- 
tion of lithium on W(112) surface, i.e., d — 1.5D, A = 
-0.137eyA, kp = OAlA-\ ip = 1.267r, Eb = -0.05eV, 
and a = 3. 16 A. We wiU use a reduced chemical poten- 
tial ^/^o with fj,o = ksT, T = lOOK. Although the 
interaction described by Eq. (|2J) is long-range we intro- 
duce finite range of this interaction in Monte Carlo study. 
The dipole-dipole interaction is accounted for distances 
between adatoms smaller than i?™ and the indirect in- 
teraction along the y-direction for distances smaller than 
(see Fig.HJ. 

The adatom coverage 6 is defined as 

i 

where N is the number of adsorption sites and angle 
brackets mean a thermodynamic average. 

It is convenient to transcript the lattice gas model as 



Ising model by introducing spin variables Si = 1 — 2ni 

id i 

where the magnetic field H is expressed by the chemical 
and interaction potentials: 

3 

Using Ising Hamiltonian, Eq. Q), we can perform sim- 
ulation in the canonical ensemble which is equivalent to 
simulation of the lattice gas, Eq. in the grand canon- 
ical ensemble. Of course, the former method is much 
easier. 



III. MONTE CARLO METHOD 

In Monte Carlo simulation we employ standard 
Metropolis algorithm to Ising Hamiltonian, Eq. on 
finite rectangular lattice of x Ly sites with periodic 
boundary conditions (PBC). Quantities such as magneti- 
sation M, specific heat Ch, and energy E are obtained 
as 

i 

E = ^<n>, (8) 

where N = LxLy is the number of spins. We will use the 
relation between coverage and the magnetisation M : 
6* = (1 — M)/2. Since we are interested in linear chain 
structures we measure also magnetisation at each row of 
the lattice 

Mj = -^^<s^u, >, (9) 

i 

where 

r|^^ = {ia.x,iay), 
Having Mj one can calculate coverage at jth lattice row 

e, = (1 - M,)/2. 

To study phase diagram we calculate the correlation 
function < SjSi > and the structure factor defined asi^ 

'^('?) = ]^ E < ^3^1 > exp(iq(rj - r;)), (10) 

where g is a wavevector from reciprocal space. This func- 
tion is used to recognize an ordered phase. 

We encountered two difficulties in examination of 
phase diagrams by Monte Carlo method: 
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FIG. 2: Ground states in the chemical potential / lattice size 
Ly plane for range of indirect interaction Ym ~ T2. Regions 
denoted by 2 and 3 correspond to phases (2 x 1) and (3 x 1). 
For detailed description of states 2a and 2b see the text. 

(1) The dependence of the results on range of the indi- 
rect interaction. 

(2) The dependence of the resulting spin structure, at 
low temperatures, on an initial spin configuration. 

In the previous papenifi we studied this model within 
the molecular field approximation (MFA) finding sev- 
eral ordered phases of {p x 1) type in phase diagrams. 
The (p X I) structure at T = 0, has every jp row (for 
j = 1,2, , Ly/p) occupied by adatoms and remain- 
ing rows are empty. Therefore, coverage in this structure 
6 = 1/p and the structure factor S{q) reaches maximum 
at q = (0 , 2tt / (pay)) . The phase diagram constructed 
for the Li/W(112) system contains the (2 x 1), (3 x 1), 
(4 X 1), (6 X 1), and (9 x 1) linear chain structures. Per- 
forming Monte Carlo simulations we observed that addi- 
tional ordered structures can occur in a phase diagram. 
It may happen when the range of indirect interaction 
Ym is shorter than the linear lattice size along the y- 
direction Ly. We start investigation of this problem by 
checking the influence of indirect interaction range on 
the ground states. As an example we present in Fig. [21 
a part of ground-state diagram for constant range of in- 
teraction Ym = 72 and different values of the lattice size 
Ly. We see that only two phases, (2 x 1) and (3 x 1), are 
present in this part of phase diagram when the interac- 
tion range Ym is close to the lattice size {Ly < 108). An 
additional phase denoted as (2a) appears for Ly > 108. 
This is linear chain structure consisting of two domains 
of the (2 X 1) phase and separated by two empty rows. 
Thus it has coverage 9 — {Ly — 2)/{2Ly). It is worth 
noting that the (2 x 1) phase has two domains: one 
with adsorbate chains placed in even lattice rows and 
the second with chains in odd rows. The configuration 
of the (2a) phase at Ly = 144 may be represented by 



sequence of distances between subsequent chains in the 
following way {235, 3i, 234, 3i}, where Ij denotes sequence 
of j identical distances I. When range of interaction 
is nearly 1/4 of Ly then another phase (26), contain- 
ing four domain walls, occurs in the phase diagram (see 
Fig. [2J|. It looks like double (2a) phase, e.g., the (26) 
phase can be represented by following sequences of dis- 
tances {234,3i,234,3i,23i,3i,227,3i} for Ly = 264. In 
a part of diagram corresponding to lower coverage ( not 
shown in Fig. ^ there are also phases containing walls 
and domains of (4 x 1), and (6 x 1) long periodic struc- 
tures but detailed discussion will not be presented here. 
It worth to mention at this place, that the(2a) phase can 
occur in the temperature phase diagram, even if the (2a) 
phase is not a ground state. In order to check the role of 
the boundary conditions, we performed additional calcu- 
lation of the ground states using free edge boundary con- 
dition along the y-direction. The same defected phases 
where obtained but the number of domain walls was re- 
duced by one. Above discussion shows that truncation of 
long-range indirect interaction in our model causes some 
technical difficulties in applying Monte Carlo method, 
e.g. it can be impossible to study finite-size scaling. On 
the other hand, the truncation of the interaction inter- 
preted as screening effect might be important from phys- 
ical point of view because it generates defected phases 
containing domain walls. 

The second problem observed in our Monte Carlo cal- 
culation, dependence of resulting configuration on the 
initial spin configuration, is due to presence of short 
range attractive interaction between the nearest neigh- 
bours along the x-direction. When linear chain structure 
is formed at low temperature, it is difficult to change 
the period of such structure by changing temperature 
or chemical potential using single spin flip Metropolis 
algorithm. We study a phase diagram performing MC 
runs at constant magnetic fleld (or chemical potential) by 
cooling down or heating up the sample. The simulation 
start from high temperature in the cooling down process 
and after thermalization and measurement temperature 
is lowering the by AT then such procedure is repeated till 
T = 0. At each temperature, all quantities are measured 
after 2 x 10"* MC steps per spin during subsequent 5000 
MC steps. The final spin configuration at T is used as 
initial one at T — AT. For some values of H the cooling 
down process finishes at T = with configuration having 
energy greater than energy of ground-state. Hence this 
way can produce metastable phases. To overcome this 
problem we use heating up process which starts at tem- 
perature AT with ground-state configuration as the ini- 
tial spin configuration. In subsequent temperature steps, 
Ti = Ti_i + AT, the specific heat Ch and the mean en- 
ergy < Ti. > are measured. This allow us to calculate the 
entropy 

S{T) = £ ^dT', 
by numerical integration. Then the free energy can be 
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obtained as 

F{T) =<n> -TS{T). 

Repeating simulations at given H with different ground- 
state configurations as initial configurations, we can find 
a phase with the lowest free energy at temperature T. 

IV. PHASE DIAGRAMS 

To study phase diagrams by Monte Carlo simulations 
we used method described in the previous section, i.e., 
simulations at constant chemical potential (or magnetic 
field) with variable temperature. The phase transition 
between ordered phases is determined by examination 
of free energies of appropriate phases. Temperature of 
transition from disordered phase to an ordered one is 
determined from behaviour of the structure factor S{q), 
specific heat Ch or sublattice magnetisation Mj. Lat- 
tice sizes Lx, Ly used in our simulations are rather small 
due to long-range interaction and large number of MC 
runs needed to construct phase diagrams. Most simula- 
tions were performed for 20 x 72 lattice sites, but some 
calculations were performed on 40 x 72 and 20 x 144. Pe- 
riodic boundary conditions were assumed, therefore Ly 
was chosen as multiplicity of 36 in order to allow forma- 
tion of ideal chain structures with period 2, 3, 4, 6, and 
9, which exist in the ground-state phase diagram. 

It has been shown in the previous section that defected 
ordered structure can occur when Ym < Ly. We now dis- 
cuss the thermal stability of such structure for two sets 
of parameters Ly and Ym- We see in Fig. |2K that the 
(2a) phase dominates the part of phase diagram corre- 
sponding to intermediate coverage and it pushes down 
the phase (2 x 1). Hence, it is impossible to pass directly 
from the ordered (2 x 1) phase to the disordered phase. 
More interesting case is shown in Fig. where the (2a) 
phase occurs at high temperatures, although it is not sta- 
ble at T=0. In this case, the (2a) phase is also placed 
between the (2 x 1) phase and the disordered phase. We 
think, that this is important result because it shows that 
analysis of the ground states is not sufficient to find all 
ordered phases. The phase (2a) has two interesting fea- 
tures 

(1) The wavevector q = (0, 27r(Lj, — 2)/ {2Lyay), at 
which the structure factor has maximum, does not 
depend on temperature. 

(2) The sublattice order parameters vanish at some 
temperatures whereas the behaviour of S{q) indi- 
cates the existence of ordered structure. 

The first feature indicates that the phase (2a) is not an in- 
commensurate phase. In order to understand the second 
feature we performed simulations with different Monte 
Carlo measurement times. If the time of measurement is 
not too long (e.g., 5x 10** MC time steps at T = 300K) we 
observe (see Fig. 0^) that order in this structure changes 
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FIG. 3: Phase diagrams in the temperature / the chemical 
potential plane for (a) Ym ~ 72, Ly = 144, (b) Ym = 36, 
Ly — 72. Monte Carlo results are represented by circles. 

spatially along y-direction. The disorder appears near 
domain walls but as one moves form the wall the chain 
structure emerges with increasing ordering. On the other 
hand we found that domain walls and the whole structure 
are moving along the y direction during simulation. This 
is a reason why the sublattice order parameters vanish if 
averaging time is long enough (see Fig. 0] where results of 
simulations with different averaging time are presented) 
but the structure factor remains finite. 



A. Application to Li/W(112) 

Linear chain structures (2x1), (3x1), and (4 x 1) 
were observed in the experiment of adsorption of lithium 
atoms on W(112) surface^^. Using MFA method, we have 
showniS. that lattice gas model, Eq. (Q), describes such 
sequence of structures. Now we present results of Monte 
Carlo simulations. To avoid occurrence of phases with 
domain walls of (2a) type, we chose the range of indirect 
interaction equal to the linear size of the lattice along 
y direction, i.e., Ym = Ly. Most simulations were per- 
formed on 20 X 72 lattice and sometimes the lattice of 
40 X 72 sites was used. The range of dipole-dipole in- 
teraction was Rm = 12ay. The ranges of these interac- 
tions used in previous MFA studieaiSi were both equal 
to 36aj^. Therefore we also repeated the MFA calcula- 
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FIG. 4: The deviation of row coverage 6j from the lattice 
coverage 9 (denoted by bars) calculated at T = SOOif, /i = 
22/io, y-m ~ 72, Ly = 144, and for averaging time (a) 5 x 10'' 
(b) 10® MC steps per spin. 




FIG. 5: Phase diagram in the temperature / the chemical 
potential plane for Ym — Ly — 72. Mean-field results are 
denoted by thick lines and Monte Carlo results are represented 
by circles (a thin line serve as guide for an eye). 



tions for present values of interaction ranges in order to 
compare the MFA results with Monte Carlo simulation 
results. As previously, MFA calculation were performed 
for all structures of {p x 1) type with p up to 12. In 
Fig. we plot results of MFA and MC methods in the 
(T,/i) phase diagram. A low temperature part of the 
diagram does not depend on the method. However in 
the upper part of the diagram situation is quite differ- 
ent. Thermal fluctuations, not present in MFA, destroy 
ordered linear chains structures at temperatures much 
smaller than MFA critical temperature (the MFA critical 




400 



300 - 



200 - 



100 - 



/ 4 




(6,4) \ 


(4,3) 


1 1 
0.2 


0.3 




e 



FIG. 6: Monte Carlo phase diagram in the temperature / 
coverage plane for Ym — Ly = 72. Results of simulations are 
represented by circles. 



temperature for the (2x1) phase can be even two times 
greater!). Moreover, MFA results do not allow to reach 
long periodic phases (4x1) and (6x1) directly from the 
disordered phase. On the other hand results of MC simu- 
lations show that all ordered phases can be reached from 
disordered phase in wide range of chemical potential. It 
is interesting to look at results of MC study presented in 
the (temperature - coverage) phase diagram (see Fig. El . 
Ordered phases of (p x 1) are separated by mixed phases 
{p,p') which denotes mixture of (px 1) phase and (p' x 1) 
phase. It is worth noting that ordered chain structures 
with period 6 and 9, present in the phase diagram, were 
not observed in experiments^ which were performed at 
temperatures above 77K. We expect that these phases 
could be observed experimentally at temperatures below 
77K. 

Now, we are going to discuss the order of the phase 
transitions. The phase transitions between the ordered 
phases, e.g., (2 x 1) ^ (3 x 1), are of the first order 
because sublattice order parameters and structure fac- 
tors change discontinuously at the transition points. The 
analysis of the phase transition from an ordered phase to 
the disordered phase is much more complicated. Check- 
ing temperature dependence of the structure factor and 
heat capacity (some results are shown in Fig.EJ) we found 
that the transition to (9 x 1) and (6x1) phases is the 
first-order phase transition because of discontinuity of the 
structure factor. The transition to the (3 x 1) and (4 x 1) 
phases seems to be continuous transition. It is difficult 
to state the order of the phase transition form the disor- 
dered phase to the (2 x 1) phase. Although, the structure 
factor seems to be continuous function of temperature, 
the heat capacity is smeared near the transition temper- 
ature indicating the first-order phase transition. On the 
oder hand, in the MFA calculation this phase transition 
is of the second order. One of possible explanation of 
this discrepancy is based on observation of additional, 
very weak maxima of the structure factor with the wave 
vector close to q = (0, 7r/(ay)) and temperature in the 
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FIG. 7: Temperature dependence of (a) the structure factor 
(b) heat capacity in (4 x 1), (3 x 1), and (2 x 1) phases for 
m/mo ~ 1-8, 6, and 19, respectively. 



narrow range of the transition temperature. This can 
suggest that defected phases hke the 2(a) phase are stih 
present close to the transition temperature. In order to 
answer whether this problem is connected with the fi- 
nite range of indirect interaction one should performed 
finite-size calculation with infinite range of the interac- 
tion. It can be done by employing a numerical method 
to accelerate the convergence of a Fourier series24. 



V. DISCUSSION AND CONCLUSIONS 

Our Monte Carlo analysis revealed that truncation of 
the indirect interaction changes a phase diagram even 
if the interaction range is very large. New phases con- 
taining domain walls occur at high temperatures between 
the disordered phase and the {p x 1) phase. This problem 
complicates study of critical behaviour via finite-size scal- 
ing. One possible solution of such problem is accounting 
for infinite range of indirect interaction. It can be car- 
ried out by employing a numerical method to accelerate 
the convergence of a Fourier series^ in the same way as 
in earlier paper^^. However, it would be time consum- 



ing way to study the critical behaviour by Monte Carlo 
simulations. On the other hand it would be interest- 
ing to extend an analysis of defected phases (e.g., (2a) 
phase) generated by truncation of indirect interaction. 
The structure factor of the (2a) phase has maximum at 
the wavevector slightly shifted from the position which 
corresponds to maximum in the ordered (2x1) phase. 
Hence it might be important in analysis of LEED pat- 
terns. 

Present MC calculations demonstrate that results of 
the mean field approximation are very bad especially at 
high temperatures. We have shown that MFA transition 
temperature from the disordered phase to the (2x1) 
phase is almost two times greater. Moreover, the high 
temperature part of the phase diagram has different 
topology than that obtained by MC simulations. How- 
ever, it is interesting to note that both method give very 
similar results at low temperatures. Present MC cal- 
culations confirm earlier obtained results™ that at low 
coverages of lithium on the W(112) the long periodic 
structures (9x1) and (6x1) are stable at temperatures 
below 77K. It would be very interesting to determine ex- 
perimentally the whole (T,0) phase diagram to check 
theoretical prediction. We are aware that present lattice 
gas model contains some phenomenological parameters, 
e.g., short range interaction Et which has influence on 
the transition temperature. We think that the better 
way to estimate the model parameters is to perform the 
first-principles calculations of the interaction potentials 
in the Li/W(112) system. 

Recently Yakovkin studied^'^ chain structures adsorbed 
on Mo(112) and Re(lOlO) using lattice gas model in 
which the indirect interaction takes nonzero value only 
for one from all lattice distances. Performing MC simu- 
lations in canonical ensemble he was able to reconstruct 
the sequence of phases (4x1) and (2 x 1) as observed 
experimentally in Li/Mo(112). It is interesting whether 
such simple model could describe correctly the phases 
observed in adsorption of lithium on W(112). 

We observed that Monte Carlo study of phase diagram 
in canonical ensemble is very difficult due to dependence 
of final MC results on initial lattice gas configuration. 
Although similar difficulty appears in grand canonical 
ensemble it is rather easy to find ground-state configu- 
rations and calculate free energies for different phases. 
However, it would be interesting to look for different al- 
gorithm than single spin flip in order to overcome the 
dependence MC results on initial configurations. 
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